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Abstract 

We present a theory for the ab-initio computation of NMR chemical shifts 
(<t) in condensed matter systems, using periodic boundary conditions. Our 
approach can be applied to periodic systems such as crystals, surfaces, or 
polymers and, with a super-cell technique, to non-periodic systems such as 
amorphous materials, liquids, or solids with defects. We have computed the 
hydrogen a for a set of free molecules, for an ionic crystal, LiH, and for a 
H-bonded crystal, HF, using density functional theory in the local density 
approximation. The results are in excellent agreement with experimental 
data. 
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Nuclear magnetic resonance (NMR) is one of the most widely used experimental tech- 
niques in structural chemistry. In particular, the chemical shift (a) spectra are a fingerprint 
of the molecular geometry and the chemical structure of the material under study. Although 
the interpretation of these spectra generally relies on empirical rules, ab-initio calculations 
of a for molecules have led in many cases to an unambiguous determination of the micro- 
scopic structure Q. So far, ab-initio calculations of a have been restricted to finite systems 
such as isolated molecules or clusters fl|]||. This is a serious limitation, because most of 
the NMR experiments are performed on liquid samples. Moreover, it is nowadays possible 
to measure a also in solids with the resolution required for structural determinations [£|||. 
E.g., a-spectra have been used for the characterization of amorphous carbon ||. In this 
letter, we present a formalism to compute, from first principles, a in extended systems using 
periodic boundary conditions. Our approach can be applied to periodic systems like crys- 
tals, surfaces, or polymers and, using a super-cell technique, to non-periodic systems such as 
amorphous materials, liquids, or solids with defects. In the case of the amorphous solid or 
liquid, the atomic coordinates may be generated by ab-initio molecular dynamic simulations 

The chemical shift measures the local magnetic field in a sample induced by a uniform 
applied magnetic field. The computation of o in an extended system is not straightforward, 
since the expectation values of the individual terms of the perturbative Hamiltonian for 
extended eigenstates are not well-defined quantities 0. To overcome this problem, we follow 
Ref. ||10|| , in which a theory for the computation of the macroscopic magnetic susceptibility 



is presented. In particular: i) we obtain the magnetic response to a uniform field as the 
long-wave limit of a periodic field; and ii) we use a generalized f-sum rule to remove the 
numerical instability which occurs in this limit. We apply the resulting equations to real 
systems, describing the electronic structure within density functional theory (DFT) in the 
local density approximation (LDA). We compute the hydrogen a for a set of small molecules, 
for an ionic crystal, and for a H-bonded crystal. Our results are in excellent agreement with 
experimental data. 



A uniform, external magnetic field B ext , applied to a sample induces an electronic current 
density Jj n (r). This current produces an induced magnetic field B in (r). If B ext is small 
enough, a condition realized in NMR experiments, then: 

B in (r) = -a(r)B ext . (1) 

Here <r(r) is the chemical shift tensor. With NMR spectroscopy, it is possible to measure 
the symmetric part of <?(r), or more often its trace, cr(r) = (l/3)Tr[er(r)], at the position of 
the non-zero spin nuclei. 

In the bulk of a periodic system, a(r) is also periodic. We may write: 

<x(r) = £?(Gy G ' r , (2) 

G 

where G are the reciprocal lattice vectors. For G^O, <5"(G) is a bulk property: 

£(G) = -47rx(G,0), (3) 

where x(G, G') is the magnetic susceptibility matrix. However, for G = 0, cx(0) is not a bulk 
property. Its value depends on the shape of the sample, and is determined by macroscopic 
magnetostatic. In our calculations we assume a spherical sample, for which: 

2(0) = -yX(0,0), (4) 



where x(0, 0) is the macroscopic susceptibility JTI ]. Thus the calculation of a in a peri- 
odic system requires the knowledge of x(G, 0). We compute the macroscopic susceptibility 
x(0, 0) following Ref. |I(J. The other elements of x(G, 0) are computed as described below. 



The susceptibility matrix is defined as the second derivative of the total energy with 
respect to the external magnetic field. In particular, 

d 2 E[B] 



b G -x(G,0)b o 



= ~K,-G, (5) 

B=0 



dB dB^ G 

where E\B] is the total energy of the system per unit volume in the external magnetic field 
B(r), B(r) = [i?obo + 5-Gb-G exp(— iG ■ r)], and b_c and b are vectors of unit length. 



Thus x(G, 0) can be evaluated using perturbation theory. However the expectation values 
of the perturbative Hamiltonians for a uniform field between extended eigenstates are ill 
defined. To avoid this problem, we modulate the external periodic field with a finite wave- 



vector q, B(r) = [£> q b q exp(zq • r) + £> g qb G^ q exp(— i(G + q) • r)] [12|], and we recover 
the results for the uniform field by considering the limit [|TTJ[] : 

b_ G -HG,0)b = -lim<_ G _ q . (6) 

We now consider a spin compensated system described by a single particle Hamiltonian. 
The derivatives of the Hamiltonian required to compute E" _ G _ q are: 

K= y (e* q r a q -V + a q - Ve iq r ) , 
#i G - q = Y (e- (G+q) r a_ G _ q • V 

+ a G q • Ve^ G+ ^ r ) , 
-^q,-G-q = a q ■ a_ G _ q e tG r , (7) 

where a s = (is x b s )/(cs 2 ), and c is the speed of light. Using perturbation theory we obtain: 

E q ,-G-q= Z( a q' a -G-q, G, q) 

+2a q ■ a G q / — — X>k,i|e~ <G '>k,i>> (8) 

where |wk,i) is the periodic part of the Bloch eigenstate of the unperturbed Hamiltonian 
with eigenvalue e k ,i, O are the sets of occupied bands, and 

/(ai, a 2 , G, q) = / d 3 r[/i(ai, a 2 , r, q) 

+r(at,a;,r,-q)]e- tG - r , (9) 



where the integral is performed in the periodic cell, with 

Ma 1 ,a 2 ,r,q) = /A^|: 

[(« kji |r)(r|a 2 -(-iV + k + q)|u^ ai ) 

+ K 4 |a 2 -HV + k)|r)(rKf)]. (10) 



I^kf 1 ) i s the first order change of the eigenstate |wk,i) due to a field with wavevector q. It 
can be obtained by solving the linear system: 

(e kji -# k+q )Kf ) = Q k+qai - f-iV + k+|) |« k)i ), (11) 

where Q k+q = (1 — J2ieo l n k+q,i)(^k+q,i|) is the projector onto the empty subspace. 

The first term on the r.h.s. of Eq. (H) is obtained as second order perturbation with the 
first order derivatives of the Hamiltonian H' and H' n „. The second term in the r.h.s. 

q Li q 

of Eq. (|8|) is obtained as first order perturbation with the second order derivative of the 
Hamiltonian i?q_ G _ q . Since a q diverges as 1/q for q —>■ 0, the two terms on the r.h.s. 
of Eq. (|§|) individually diverge as 1/q. To remove this divergence, which would produce a 
numerical instability, we use the generalized f-sum rule: 

7^3 EK*l e >M> = -/(ai,a 2 ,G,0). (12) 

Substituting the f-sum rule into Eq. (§) we obtain: E" G = [/(a q , a_G- q , G, q) — 
/(aq, a_ G _ q , G, 0)]. Then, for G ^ 0, 



b_ G -x(G,0)b = -lim^_ G _ q 
d 

-f(q x b ,G x b_ G ,G,gq) 



(13) 

g=0 



c 2 G 2 dq 

where q is the unit vector in the direction of q. Finally, the derivative with respect to q in 
Eq. ( |13D can also be evaluated using the following limit: 

b- G • X(G, 0)b = - lim[/(q x b , G x b_ G , G, qq) 

-f(q x b , G x b_ G , G, -qq)}/(2qc 2 G 2 ). (14) 

Note that, for G^O, x(G, 0) is proportional to the first derivative of / with respect to q, 
whereas the macroscopic susceptibility %(0, 0) is proportional to the second derivative of / 
with respect to q fllPf . 

In practice, we evaluate numerically x(G, 0) using Eq. ([H]) with a small, but finite q, and 
the k-integral in Eq. (|It]) with a finite summation in the irreducible wedge of the Brillouin 
zone. 



We describe the electrons using DFT-LDA; i.e., we neglect any explicit dependence of 
the exchange-correlation functional (E xc ) on the current density. The current dependence 
of E xc could be taken into account using the approximate functional proposed in Ref. flT5 



but in practice, this produces only negligible corrections to a in real systems ||. An ad-hoc 
procedure to include many-body effects beyond DFT in the the calculation of a has been 



proposed in ||14|| . While this approach improves over DFT in small molecules, the corrections 
to DFT vanish for periodic systems, where the eigenstates are always extended. In general, 
to compute the second order variation in the DFT total energy with respect to an external 
perturbation, one should take into account the linear variation of the Hamiltonian induced 
by the linear variation of the charge bp. However, if the perturbation is a magnetic field, bp 
is zero by time reversal symmetry. Thus Eq.s (0-0) are correct within DFT. 

In the present calculation, we will consider the magnetic response of valence electrons 



only. We describe the ionic cores by norm conserving pseudopotentials |15] in the Kleinman 



Bylander form This approximation does not affect a of the nuclei without core electrons, 
such as H PUTT], but those containing core electrons. In the latter case a computed with 
pseudopotentials differs from the one computed with an all-electron scheme by three dif- 
ferent terms: i) the diamagnetic core contribution, which is independent of the chemical 



environment; ii) a contribution due to the transitions from valence states to core states |L0 
iii) a contribution due to the difference between the all-electron valence wavefunctions and 
the pseudo wavefunctions in the core region. We found that, for first row atoms like carbon, 
the error due to the pseudopotential is minor, since the terms ii) and iii) are usually much 
smaller than the range of variation of a with the chemical environment. In the present 
paper, however, we present only the results for a of H, which is not affected by the use of 
pseudopotentials. Finally, in our pseudopotential calculation, we replace (— iV + k + q/2) 
in Eq. (|TT| ) with (v£ + v£ +q )/2, where v£ = (d/dk)H£ is the velocity operator, and if£ is 
the pseudo-Hamiltonian [|18j. 

We expand the wave-functions using a plane wave (PW) basis set. We obtain Iw^f 1 ) 
by solving Eq. ( |Tl~l) with a conjugate gradient minimization |19 |; i.e., we never use the 
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unoccupied eigenstates of H^. We obtain /(ai, a 2 , G, q) from Iw^T 1 ) using fast Fourier 
transforms. Thus, the computational effort to evaluate a for all nuclei in the sample is 
comparable to that required for a computation of the total energy P0| . We compute o in 
systems containing atoms of H, C, Li, and F. For H we use a local pseudopotential to speed 
up the convergence with respect to the PW basis PTJ] . For C, Li and F atoms we use a 



pseudopotential with a non-local s-projector, to describe the interaction of valence states 
with the Is core states. We use a PW energy cut-off of 100 Ry for the F atom, and of 70 Ry 
for the other atoms. Molecular and crystalline structures are taken from experiments |f22| . 



The results for the free molecules are obtained using super-cells, such that the error on o 
due to interaction between a molecule and its periodic replica is smaller than 0.1 parts per 
million (ppm). 

The results for a set of molecules are reported in Table |I[ All the computed values are 
in excellent agreement with experimental data. The discrepancies are larger for C2H4 and 
C2H2, in which a double and a triple C-C bond exist, respectively. This could indicate a 
relatively larger inaccuracy of LDA in describing these types of bonds. Similar agreement 
with experiments for the hydrogen a in molecules is found in previous calculations 

The results for molecular and crystalline LiH are reported in Table y. LiH is an ionic 
system, which crystallizes in the rock-salt structure. The Li-H equilibrium distance in the 
crystal is significantly larger than that of the molecule. Our results show that the difference 
between a in the molecule and in the solid is very small. To verify if a is sensitive to the 
geometry, we consider a crystal at a pressure of 65 GPa, in which the Li-H distance is equal 
to that of the molecule. In this case we find a much larger value of a. Thus, a of the crystal 
and of the free molecule are only similar at the equilibrium geometries. 

The results for molecular and crystalline HF are reported in Table [TIT]. The intra- 
molecular bond of HF is covalent. In the liquid and solid phases, the HF molecules are 
bonded together in zigzag chains via H-bonds. HF has the strongest H-bonding found in 
nature. This strong inter-molecular bond is reflected in the large variation of a observed in 
the transition from the gas phase to the liquid phase, and in the large temperature depen- 



dence of a in the liquid phase |24| . The HF crystal is constituted of a 2D stacking of parallel 
ID zigzag chains. We compute a for the crystal at the experimental geometries measured 
for a DF crystal at two temperatures, 4.2K and 85K [P5fl . The computed a in the free HF 
molecule is in outstanding agreement with experimental data. In the crystal, the computed 
a decreases dramatically to a value close to the one measured in the liquid phase (we are not 
aware of any experimental measurement of a for solid HF). Also the theoretical variation of 
o with temperature, 8a /ST = 0.013 ppm/K, is very close to the one measured in the liquid, 
which is da/dT = 0.0135 ppm/K f24|. Thus, with respect to a, the behavior of the solid is 
very similar to that of the liquid. Interestingly, a single H-bond does not account for the 
large decrease of a observed in the condensed phases. Indeed the computed a for the central 
H atom in a HF-HF dimer, is closer to the a of the free molecule than to that of the liquid 
i- 

In conclusion, we have presented an ab-initio theory for the evaluation of NMR chemical 
shifts (cr) in extended systems. We have computed, for the first time, a in real solids. Our 
results show that DFT-LDA predicts hydrogen a which are in excellent agreement with 
experiment. The evaluation of a requires the same numerical effort as the computation of 
the total energy PU| . Thus, our approach can be applied to super-cells containing up to a 



hundred atoms, which are sufficient to model liquids or amorphous materials 
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TABLES 



cr-theory (ppm) u-experiment (ppm) 



H 2 25.9 26. 2 a 

CH 4 30.7 30.61 a 

C 2 H 6 29.7 29.86 fe 

C 2 H 4 24.5 25.43 6 

C 2 H 2 28.6 29.26 6 

TABLE I. Computed and measured H chemical shifts (a) for a set of free molecules. Ref. 
[23|, b Ref. 



d-theory (ppm) 

LiH molecule do = 1-595 26.6 

LiH crystal a = a = 4.08 26.3 

LiH crystal a = 2d = 3.19 31.2 

TABLE II. Computed H chemical shifts (a) in LiH for a free molecule and a spherical crys- 
talline sample (in the rock-salt structure), do is the equilibrium bond- length of the molecule, a 
is the lattice constant, and ao is the experimental equilibrium lattice constant. For a = 2do, the 
shortest Li-H distance in the crystal is equal to do. 
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T (K) d-theory (ppm) c-experiment (ppm) 

HF molecule 28.4 28.5 a 

HF crystal 4.2 20.71 

HF crystal 85 21.76 

HF liquid 214 21.45 6 

TABLE III. Computed and measured hydrogen chemical shifts a in HF for a free molecule, a 
spherical crystalline sample, and a spherical liquid sample. T is the temperature. a Ref. Ref. 



13 



